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ABSTRACT 

The methods of obtaining the average spectral shape in a low statistics regime are 
presented. Different approaches to averaging are extensively tested with simulated 
spectra, based on the ASCA responses. The issue of binning up the spectrum before 
fitting is discussed together with the choice of statistic used to model the spectral 
shape. The best results are obtained with methods in which input data are represented 
by probability density functions. Application of weights, representing the coverage 
between the input and output bins, slightly improves the resolution of averaging. 

Key words: spectral shape - Fe Ka line. 



1 INTRODUCTION 

Spectral information obtained from X-ray astronomy instru- 
ments is usually the result of a compromise between the 
aim of achieving the best possible energy/spatial resolution 
and the purpose of collecting a large number of photons. 

, In consequence, there will be always some class of objects 
which are too faint, with spectra which cannot be studied 
with all desired particulars. The problem of a lack of statis- 

, tics can be solved to some extent by averaging a number of 
weak spectra, but this implies examination of rather com- 

\ mon properties. Nevertheless, since spectra can be grouped 
into subsamples according to some better established quan- 
tities such as the continuum slope, flux, hardness ratio, etc., 
this method may be quite powerful for studying various cor- 
relations. 

A clear distinction should be made between the aver- 
age spectrum and the average spectral shape. The former is 
simply the average flux and the result is dominated by the 
brightest objects or states (in the case of average for a single 
object). The latter is the average of a relative quantity, spec- 
tral shape, usually defined as the ratio between data and a 
simple continuum model, common for all studied objects. 
Such a proportion is often employed to bring some discrete 
features into prominence — a well known example is fig. 1 
of Tanaka et al. (1995) where the redshifted iron Ka line 
profile, observed for MCG-6-30-15, was shown. The shape 
defined as above is customarily used only for illustrative pur- 
poses, e.g., Reynolds (1997) and Reeves (2003). However, the 
average shape can be constructed and studied in a quantita- 
tive way, as was done for the spectra of Syl nuclei observed 
by ASCA (Nandra et al. 1997a,b). A similar investigation 
was performed later for a larger sample of Syl ASCA ob- 
servations (Lubihski & Zdziarski 2001), where average shape 
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spectra were obtained for subsamples grouped according to 
the continuum slope. The 'average shape spectrum' is de- 
fined as the average data to continuum model ratio (i.e., the 
average shape), multiplied by the average continuum model 
(model with the average parameters calculated with a stan- 
dard weighted mean). 

The average profiles of the iron Ka line presented in 
Nandra et al. (1997a) and in Lubinski & Zdziarski (2001) are 
clearly distinct. This difference was ascribed partially to the 
change in calibration of ASCA SIS detectors done after 1997 
and partially to a different approach to averaging. However, 
the first explanation was later questioned and the whole dif- 
ference was assigned to the dissimilar averaging procedure 
(Yaqoob et al. 2002). The issue of changed calibration will be 
discussed elsewhere; here we want to consider the problem 
of correct averaging. Additional motivation comes from the 
fact that it seems promising to apply similar averaging pro- 
cedures to the data from other missions, such as Chandra, 
X M M—N ewton, and, in future, Astro-E. Therefore, it is 
important to have at one's disposal a verified method, ex- 
ploiting all available information in the most efficient and 
accurate way. 

In the following sections we discuss three basic aspects 
of spectral shape averaging: bin weights, prebinning and the 
character of the data. The first is the way in which the infor- 
mation on the relative positions of the input and output bins 
is taken into account. Prebinning means here the summing of 
counts from a single input spectrum over the span of the out- 
put bin. Finally, the Poisson character of the data is impor- 
tant for low numbers of counts where the standard weighted 
average, which assumes a symmetric, Gaussian probability 
density distribution, cannot be used. These aspects are con- 
nected — for example, prebinning leads to the loss of some 
information on the data distribution in input bins, but, on 
the other hand, increases the number of counts. Various ap- 
proaches to these basic issues can be combined to construct 
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different averaging methods; our aim is to test them through 
extensive tests performed on simulated data. 



2 AVERAGING METHOD 
2.1 Rebinning 

By definition, the average value of a function £(e) over the 
range of its argument, (e», e»+i), is equal to the ratio of this 
function integral over this range to the length of this range, 
Ae,, 



f:; +1 q(e)de _ ge)de 



(1) 



Here and after we use the following convention: angle 
brackets (} denote the average over some range of the func- 
tion argument, i.e., rebinned value, without weights associ- 
ated with the accuracy of the averaged quantity. Any aver- 
age, for a single spectrum or for many spectra, weighted by 
the accuracy weights is indicated by a dash over the symbol. 
Subscripts i,j,k are used for input data, spectra and out- 
put data, respectively. To distinguish if the input data are 
summed or averaged for a single spectrum or for all spectra, 
we will use different upper summation limits, n k j for a sin- 
gle, j-th spectrum and n k s for all spectra. The output data 
for an individual spectrum will be denoted by an additional 
index j. Finally, the number of averaged spectra is equal to 
n s . 

Assume that we know the averages of a certain energy 
function £ for some initial distribution of energy ranges Aa 
and we want to determine averages for another set of en- 
ergy ranges AE k - This is rebinning: for a given spectrum, 
values in some bins are converted to values in bins occu- 
pying different ranges. The idea of rebinning is illustrated 
in Fig. 1. In general, the output bin k, with boundaries 
(Ek, .Efc+i), for a spectrum numbered with j expands over 



rikj input bins (with lower limits eo, ei, 



and upper 



limits ei,e2, ...,e„ fe , ). The first and the last input bins can 
lie partially outside the output bin. Input data are discrete, 
we know only the average of the unknown function £(e) mea- 
sured by the detector over the i-th bin, its value is at- 
tributed to the centre of the bin. Using the approximation 
that £(e) is constant within an input bin and introducing 
bin weights bi we obtain the formula for the rebinned value 
{Qkj in the form 



<0w = J>® i - 

where 

gi - E k 
AE k ' 
Aej 
b, = { A~EJ' 

Ek+i — e n . 



(2) 
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partial overlap, left boundary; 
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partial overlap, right boundary; 
input bin covers output bin, 



(3) 



and 



5> = i. 



X 




Figure 1. Idea of rebinning of function over some range of its 
argument. The average for segment (E k , E k+ i) is determined 
using the known averages for segments (eo,ei), (ei,e2), 



-l,en 



(4) 



Defining 5, as the width of the overlap between input 
bin i and output bin k, the weights bi may be simply ex- 
pressed as 5i/AE k . 

Boundary input bins, with the same value of weight bi, 
can occupy quite different ranges outside the output bin. 
Hence, they can represent different information on the aver- 
aged function, integrated over a different range of argument. 
To take this fact into account one can apply another bin 
weight, inversely proportional to the input bin width, equal 
to 1/Ae» or, better, to 6»/Ae». 

Figure 2 shows an example of averaging with and with- 
out bin weights. The averaged function f(E) is a constant 
plus Gaussian peak plus an edge, modelled by a negative 
half- Gaussian. The widths of Gaussians were set to 0.02 
keV and 0.1 keV for peak and edge, respectively. This func- 
tion was integrated over 0.01 keV bins to obtain the refer- 
ence shape after using the instrument with 0.01 keV energy 
resolution. The function was then integrated 200 times for 
random binning patterns, with bin widths Ae^ taken from 
some interval. The resulting 'spectra' were rebinned in four 
different ways (with bin weights listed in Fig. 2) to 0.01 
keV output bins and the results were compared with the 
reference shape. The whole procedure was repeated for dif- 
ferent ranges of input bin widths. When the input bins Aa 
are much broader than 0.01 keV the resulting averages with 
and without bin weights are indistinguishable. For input bin 
widths comparable to or less than 0.01 keV, the differences 
in results correlate with decreasing input bin width. Proce- 
dures with bin weights clearly better reproduce the sharp 
features of the averaged function than the simple arithmetic 
average, the best of these is the weight in the form &f/Ae» 
but the differences are rather small. 

Bin weights play a role similar to approximating the 
shape of the rebinned function with a polynomial spline. 
However, for poor quality data they are safer than any iter- 
ative spline procedure, because the latter method may am- 
plify some spurious spectral features during consecutive iter- 
ations. Therefore, since we have to work with low statistics 
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E [keV] 

Figure 2. Upper panel: Averaging of function with and without 
bin weights. Solid line shows the averaged function f(E), reference 
shape (dots) is the result of integrating this function over 0.01 keV 
bins. Averaged data are obtained by integrating f(E) with random 
bin widths taken from 0.01-0.02 keV range (an example of these 
bins is shown with the horizontal lines). The average shapes were 
calculated for 200 input data sets, with output bin widths equal 
to 0.01 keV and using 4 different bin weights. Lower panel shows 
the differences between the results of averaging and the reference 
shape. 



data and since the energy resolution of ASCA SIS detec- 
tors is not so high, the discrete spectral features can be well 
traced by bin weights alone, without using an additional 
spline approximation. Moreover, any more complex spline 
procedure applied to a small set of weak spectra may lead to 
a quite accidental approximation of the local spectral shape. 
On the contrary, bin weight values are well defined, do not 
depend on the averaged function value and control the rel- 
evance of information given by the input bin, comparing 
simply input and output bin locations. 



2.2 Prebinning 

The quantity which we intend to average is the ratio of mea- 
sured and modelled fluxes. Since the flux is proportional to 
the number of counts we can utilize the fact that the num- 
ber of counts measured for a broader bin is equal to the sum 
of counts collected for narrower bins constituting this broad 
bin. Then, instead of averaging input flux ratios for a given 
output bin, we can directly determine data and model fluxes 



for that bin. The measured flux, (f)kj, is equal to the sum of 
all net counts D t , divided by appropriate detector area Ai, 
and normalized with the observation time, T, and output 
bin width, AE k , 



if) 



kj 



TAE k 



(5) 



where weights gi — Si/Aei (= biAEk / Aa) are intro- 
duced for boundary input bins, only partially overlapping 
with the output bin. 

In practice, there is no simple method of performing 
spectral fitting with arbitrary energy bins, i.e., with energy 
bins adjusted to cover the output bins. This is due to the 
fact that the instrumental response function is defined as a 
matrix for a fixed set of energy bins. Therefore, it is neces- 
sary to perform prebinning with model counts determined 
for input channels instead of fitting them for output bins. 
The modelled flux, (/}J™, is expressed through the model 
net counts, Mi, 



(f)Tj = 



TAEk 



(6) 



The numerator in (5) can be replaced by ^™=i /Ak, 
where the area A k represents the effective detector efficiency 
for bin k. After similar replacement in (6), the ratio of the 
data and model fluxes for output bin k is equal to 



(r)kj = 



(7) 



The procedure with prebinning is simpler than that 
based on averaging all input ratios (r), (=Dj/M»); after pre- 
binning one has only to average ratios {r)kj obtained for the 
output bins from various spectra. Nevertheless, there is one 
disadvantage: the information included in the input ratios 
{r)i is lost. The issue of how this affects the results by re- 
ducing the resolution will be discussed later, in Sec. 5. 



2.3 Accuracy weights 

The initial ratios (r); are measured with some finite accuracy 
and this should be taken into account in averaging. For the 
set of n ks independent ratios, described by the probability 
density functions Pi(r), the mean value (centre of gravity) 
is given by the integral 



r k = 



J^rpki^dr 
J^Pk(r)dr ' 



(8) 



where the joint density function Pk{r) is equal to the 
product of partial densities 



Pk(r) = JJpi(f). 



(9) 



The standard deviation for rk, <J- k , is calculated from 
the variance definition 



JZoO - r k ) 2 Pk(r)dr 

/Ho p* (»■)*■ 



(10) 
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2.4 Standard average 

In the case where the densities Pi(r) have Gaussian shape, 
centered at {r)i and with width parameters A(r)i, one 
obtains from (8) and (10) the standard formula for the 
weighted mean and its uncertainty 



Th = 



ELI w *( r )* 



Ar fe 



VET 



(11) 



where weights Wi are equal to l/(A(r)i) 2 . 

The above formula can be applied to ratios obtained 
from larger numbers of counts, when the Poisson probabil- 
ity function associated with the data can be approximated 
by a Gaussian function. In such a case, the unknown uncer- 
tainty of the true number of counts can be approximated by 
the square root of the measured number of counts. Then, 
the ratio uncertainty A{r)i is equal to \fgiNi + gtBi/giMi, 
where iV; denotes the source (net effect + background) num- 
ber of counts and Bi is the background number of counts. 
Weights gi modify these numbers for an input bin lying on 
the boundary of the output bin. In Appendix A it is shown 
that for relative quantities such as flux, the accuracy weights 
Wi to some extent play a role similar to the bin weights. 

Obviously, the condition of large number of counts can 
be more easily fulfilled in the case of prebinned data, when 
the output bins are broader than the input ones. The uncer- 
tainty of ratio {r) k j, given by equation (7), is then approxi- 
mated by 



A(r) fc3 



(12) 



The final ratio r k and its uncertainty Ar k for prebinned 
data, averaged over all spectra, is calculated from (11), with 
(r)i replaced by {r) k j, weights w k j equal to l/(A(r) fe j) 2 and 
summation going from j = 1 to j = n s . 

2.5 Combined weights 

According to the results of Sec. 2.1, the accuracy weighted 
average for input ratios {r)i should be modified to incor- 
porate bin weights. Application of another type of weight, 
such as bin weights bi, can be realized via broadening of the 
probability density distributions associated with the data 
by raising them to a power equal (or proportional) to this 
additional weight. It should be stressed that this procedure 
is used only to change the relative widths of these distri- 
butions, by taking into account the overlap between given 
input bin i and output bin k. There is, then, some arbitrari- 
ness in defining the broadening power index. In tests with 
simulated spectra, we found that the change of bi to 6i/10 
or 10&i affects only the fifth digit in the result of averag- 
ing. Nevertheless, any broadening of the initial distributions 
leads to a change of the width of resulting distribution, hence 
affects the error of r k . For narrow input bins (bi <g; 1) the fi- 
nal distribution can be quite broad and its dispersion should 
be renormalized, since bin weights are used only to redefine 
the calculation of the average, and should not change the 
accuracy of the input {r}i values. To avoid absolute narrow- 
ing of any of initial distributions the renormalization is done 
by replacing bi in the above equation by b\ — bi/max(bi). 
This leads to only slightly broader distributions than those 
obtained without bin weights. 



The combined weights for the standard average are 
equal to biWi, since the width of broadened Gaussian is equal 
to Oil \fbi. However, the combined weighted average and its 
propagated error, calculated from (11), can be expressed us- 
ing weights bi directly, due to the fact that the normalizing 
factors l/max(bi) in numerator and denominator cancel out 



Tk = 



£?=; b i W i 



Ar k = 



E"=i biWi 



(13) 



The above formulae were used to obtain the average 
spectral shape for Syl active nuclei observed by ASCA 
(Lubinski & Zdziarski 2001). However, some correction 
should be made: weights Wi for border input bins should 
be decreased by a factor equal to weight gi, since the (r)i 
error is increased by \f\Jg~i due to only a partial share of 
input bin flux in the flux of the output bin. In consequence, 
we obtain 



r k = 



Ar ks = 



(14) 



Er=i hi 9i w * ' £r=i bi 9 iWi 

As discussed in Sec. 5.2, the above correction affects 
mainly the ratio errors, leaving the ratio values almost un- 
changed. 

2.6 Average for Poisson data 

The formula given by (11) is derived from the definition of 
the mean (Eq. 8) for quantities described by the probability 
density function in Gaussian form. However, the numbers 
of counts collected in a single channel by SIS instruments 
of ASCA for weak sources such as AGNs are very small 
- above 7 keV these numbers are often equal to 1 or 0. 
Therefore, equation (11) cannot be used for a low number 
of counts, where the Poisson distribution differs substan- 
tially from the Gaussian one. Moreover, in this situation the 
unknown uncertainty of the true number of counts cannot 
be approximated by the square root of the measured number 
of counts. 

Individual densities in the counts space, jJi(A), includ- 
ing the case of border input bins with numbers of counts 
modified by weights gi, can be expressed as Poisson func- 
tions of unknown mean A, with observed giNi source counts 
and giBi background counts as parameters 



Pi (\) = Ci- 



-(A+ 9i B,) 



{X + giBi) 



9iNi 



V{g l N l + 1) 



(15) 



where constant d ensures proper normalization of the 
probability function. The above formula has the form of a 
Poisson distribution, however, here the problem is inverted 
— we are interested in a function of continuous argument A 
for given numbers of observed counts. 

Because numbers of counts modified by weights gi are 
not necessarily the integer numbers, the standard factorial 
present in the Poisson distribution definition, in Eq. (15) 
is replaced by the complete gamma function F(giNi + 1). 
For the same reason we cannot give an analytical expression 
for the normalizing factor, (cf. Eq. (B3) in Appendix B), 
because for non-integer giNi the expansion of (A + giBi) 9iNi 
binomial is not finite, and value of d has to be computed 
by numerical integration. 
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The formula (15) is based on the assumption that the 
background number of counts is known precisely. Evidently, 
in a real situation we do not known the true background 
rate. However, as will be shown in Sec. 5, this assumption is 
sufficient to give proper averaging results. Our main goal is 
to describe the data probability distribution in a way more 
valid than a Gaussian distribution, and Eq. (15) leads to a 
satisfactory result. 

Nevertheless, in Appendix B we present two solutions 
to the problem of the unknown background which we have 
found in literature. As can be expected, these procedures 
lead to a more diffuse distribution than that given by (15), 
especially with a broadened left tail when the observed back- 
ground number of counts is close to zero. We have tried to 
include these methods in the averaging code, but, because 
of some problems with the implementation of these meth- 
ods in current versions of popular spectral fitting codes (see 
next Section) we cannot model the continuum in a corre- 
sponding way. Therefore, the treatment of data in averaging 
and in continuum modelling cannot be consistent, and, in 
consequence, the data/model ratios will diverge. 

Direct calculations based on Eq. (15) are impractical 
for larger values of giNi. Therefore the function pi(X) is 
calculated using the algorithm invented for computing the 
binomial distribution (Loader 2000). Computation, except 
for the trivial case when X + giBi = 1, reduces to calculating 
the exponent of some function ai(\), depending on gi and 
Bi, normalized with some factor j3i depending on Ni, Bi and 
gi. Then the probability density can be given in the form 
e «i(A) 

P ' (A) = WW ' (16) 

Turning to ratios, we take into account the modelled 
number of counts, giMi, and the averaging reduces to deter- 
mining the mean of the joint density distribution function 



Pk(r) = Y[C ie a 



(17) 



where the function a'i is equal to a, with an argument 
scaled by the factor l/^Mj and the normalizing factors C; 
are functions of Ni, Bi, Mi and gi. 

In the case of prebinned data, for a spectrum numbered 
with j we have function a' k j and factorial C k j depending on 
the summed counts Y^i=i 9* N i, Yl™=i 9i B i and Ym=i 9i M i, 
and the joint density distribution function is calculated for 
n s spectra 




Figure 3. Four probability functions (Eq. (15)) for net counts 
D = 1 and background counts B equal to 0,1,2,3 (from top). 
For N = 1 the example of broadening of the probability density 
function by the bin weight = 0.75 is shown with a dashed line. 




Figure 4. Example of the joint probability distribution applied 
to determine the mean ratio. Thin lines show the Poisson proba- 
bility functions representing the averaged input ratios, thick solid 
line shows their product distribution, normalized to have maxi- 
mum equal to 1. Dashed line illustrates the product distribution 
obtained when the input data are given by the Gaussian proba- 
bility functions. 



(18) 



At last, when bin weights bi are taken into account for 
data without prebinning, the formula for the joint probabil- 
ity density functions has the form 



Pkif) 



n ks 



(19) 



There is no computer on Earth able to directly multiply 
many very small numbers, which are unavoidable for a broad 
range of ratio values, especially when n ks is of the order of 
thousand and the multiplied distributions are narrow. Thus 
it is necessary to replace Pk{r) with its logarithm 



log Pfc (r) = logJj[c i e a i (r) ] 6< =2log[c i e a J 

i=l »=1 

= ^bta'iir) & log Ci. 



M 



(20) 



The second sum in the last line of the above equation 
can be dropped since it only corresponds to adding a con- 
stant to the density distribution function and the computed 
density function is 



logp fc (r) = y^6j« ■(?-). 



(21) 



The final average r k is found from equation (8) as the 
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mean of the exponent of the above function but with the 
lower integral limit equal to 0, since the number of counts 
(or ratio) cannot be negative. Accordingly, the accuracy of 
r k is calculated from Eq. (10). 

The joint density function Pk(r) is constructed as a 
likelihood function in the maximum likelihood estimation 
method. Then, the best estimator for Tk is the mode, for 
which the accuracy should be calculated using the second 
derivative of Pk(r). The computation of mode for density 
function given by Eq. (8) is simple but numerical determi- 
nation of its second derivative and then its expectation value 
to assign the accuracy to r k for many input bins may signif- 
icantly increase the computation time. Therefore, because 
the final distribution functions are usually quite symmetric, 
especially for a larger number of averaged spectra, the mode 
and its error were replaced here by the mean and standard 
deviation. As it was tested, even for an extreme case of a sin- 
gle, weak spectrum, the difference between mode and mean 
of Pk{r) does not exceed 0.5%. 



3 CONTINUUM MODEL 

Considering the issue of averaging the data/model ratios 
in Sec. 2.2, we have concluded that the continuum model 
should be fitted to the non-binned spectrum, i.e., collected 
with original, single SIS channels. Due to this fact, for a 
low number of counts we also have to treat the spectral 
modelling in a non-standard manner. Since the x 2 -statistic 
cannot be applied during fitting, an alternative approach is 
needed. A fairly popular and well established solution is the 
C-statistic (Cash 1979). It is based on the maximum likeli- 
hood method, and has no limitation according to the num- 
ber of counts. In the case of Ni source counts and Si model 
counts it is defined as (Freeman, Doe & Siemiginowska 2001) 

C = 2^2(Si - NilnSi), (22) 
i=i 

where n c is the number of data channels. 

The difference of two Poisson distributions is not Pois- 
son distributed, thus the C-statistic cannot be applied to 
background subtracted data. Therefore in the C-statistic 
case the source and background spectra were fitted simulta- 
neously. The net number of model counts Mi was determined 
by subtracting the appropriate numbers modelled for a given 
channel. In consequence, data in modelling were treated in 
a way consistent with averaging based on equation (15). 

The C-statistic is implemented in both most popular 
X-ray spectra fitting codes, XSPEC (Arnaud 1996) and 
SHERPA (Freeman, Doe & Siemiginowska 2001). However, 
used within XSPEC it produced biased results: the fitted 
power law indices and normalizations were much larger than 
those assumed in the simulated model. 1 Therefore, in the 
case of non-rebinned spectra, fitting the model of continuum 
used to normalize observed data to obtain the spectral shape 
was done only with SHERPA. For comparison, we have also 

1 There was a bug in XSPEC, repaired in the its version 11.2.0bs, 
after preparation of this paper. Now XSPEC used with the C- 
statistic produces results consistent with those obtained with 
SHERPA. 



binned up spectra and fitted them with a x -statistic using 
both XSPEC and SHERPA. All tests were performed with 
XSPEC, version 11.1 and SHERPA, version 2.3. 

There is a Bayes statistic option in SHERPA, which 
corresponds to the method derived by Loredo (1992). We 
have tested also this approach during spectral fitting, un- 
fortunately the results were similar to the results obtained 
with the C-statistic and XSPEC: the fitted power law was 
much steeper than the model assumed in the simulation of 
the spectrum. Hence, since the usage of the Bayes statis- 
tic in SHERPA seems to be somewhat uncertain and since 
the results obtained with C-statistic and SHERPA are quite 
satisfying, we do not use the Bayes approach in the tests 
described later. 

All continuum models for the source and background 
spectra were fitted in the energy range (3-4.5,7.5-10) keV, 
i.e., the reference continuum shape was determined in the 
vicinity of the Fe K a line but without the line region itself. 
We used a power law model, assuming independent slope 
parameters for source and background and independent nor- 
malization for each spectrum. 



4 TESTS 

4.1 Reference shape 

The procedure developed for averaging the spectral shapes 
has to be verified, moreover, some tests of its alternatives are 
needed. It is obvious that such tests cannot be done for real 
data, as we must know the actual average shape to have 
a reference. For this purpose we have simulated 100 spec- 
tra for each of two SIS instruments. The reference model 
was similar to that obtained for the average ASCA Syl nu- 
clei spectral shape (Lubinski & Zdziarski 2001). Simulations 
were done with responses from different periods of ASCA 
mission. The width of a single SIS channel was equal to 14.6 
eV, since we used the BRIGHT2 data mode with 1024 chan- 
nels. As the starting conditions we adopted exposure times 
and background spectra obtained for five observations of 
IC 4239A, which were short, with an elapsed time of 7-18 
ks. The background spectra were extracted for rather small 
regions with a radius of 28.5 SIS pixels. Hence, we tested 
a rather extreme case of low statistics to be sure that the 
method behaves well on the boundary of its application area. 

In the upper part of Figure 5 we have shown our ref- 
erence model, consisting of a broad, disc line component, 
a narrow, Gaussian line component and continuum in the 
form of power law. The parameters of the model were as 
follows: power law index, T = 1.8, power law normalization, 
A = 2.333 x 10~ 2 keV _1 cm _2 s _1 , disc line energy, 6.4 keV, 
inner disc radius, 6 GM/c 2 , outer disc radius, 1000 GM/c 2 , 
disc emissivity in the form (1— ^/6/R)/R 3 , inclination, 45°, 
Gaussian line energy, 6.4 keV, Gaussian line width, 0.01 keV. 
The normalization of disc line and Gaussian components was 
adjusted to get equivalent widths equal to 130 eV and 60 eV, 
respectively. 

Spectral shape is defined as the ratio between observed 
data and fitted continuum model, thus the reference model 
should be transformed in the same way. The reference data 
were simulated with a very long exposure time, 10 9 s, using 
responses of original IC 4329A spectra. Then, the contin- 
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Figure 5. (a) Model function (solid line) consisting of a disc line 
and a Gaussian together with the power law continuum. Dashed 
line shows the continuum model, i.e., power law fitted outside the 
Fe line region, in the energy ranges 3-4.5 and 7.5-10 keV. (b) Ratio 
of the line and continuum model functions (solid line) compared 
with the data/model ratios (dots) obtained for the simulated ref- 
erence spectra. 



uum for the reference spectra was modelled with a power 
law model outside the iron line region. Resulting data/model 
ratios were averaged with weights equal to the original obser- 
vation times, to take into account changes of SIS responses 
during satellite operation. These reference ratios are pre- 
sented in the lower part of Fig. 5. Rebinning with finite 
bin width always leads to a change in the function shape, 
thus the results of averaging done with tested methods are 
compared not only with the reference ratios obtained for a 
single input bins, but also with the results of averaging these 
ratios to appropriate output bins. This reference averaging 
was done with our best method la (described below), but for 
such a high number of counts all procedures produce almost 
indistinguishable results. 



4.2 Tested methods 

Various averaging procedures can be constructed on the ba- 
sis of three elements: accuracy weights, bin weights and pre- 
binning. We have tested some of these combinations to check 
their behaviour in different situations. Table 1 presents the 
characteristics of the tested procedures. Methods numbered 
with I are based on determination of the joint density func- 
tion, methods II are those using the accuracy weights u>i, 
and method III is the arithmetic average with prebinning. 
Among them, method He uses the standard weighted aver- 
age (11) applied to only those input bins, whose centres lie 
in the given output bin. Since this procedure is probably the 



Table 1. Averaging methods tested with simulated spectra. PF 
denotes the probability density function. Bin weights equal to 1 
mean that the input bin is taken into account only in calculating 
the average for the output bin containing its centre. 



Method 


Accuracy 


Bin 


Prebinning 


Formula 




weights 


weights 






la 


PF 


9i 


yes 


(18) 


lb 


PF 


9i 


no 


(17) 


Ic 


PF 


bigi 


no 


(19) 


Ha 


Wi 


9% 


yes 


(7,11,12) 


lib 


Wi 


9i 


no 


(11) 


lie 


Wi 


kgi 


no 


(14) 


lid 


Wi 


bi 


no 


(13) 


He 


Wi 


1 


no 


(11) 


III 


1 


9i 


yes 


(2,7) 



one most commonly used, it will be termed 'standard' in the 
rest of the paper. 

Tested methods were applied to three types of 
'data/model' data. The first one was obtained for spectra 
with original SIS channels and the continuum model fitted 
with the C-statistic. The data of the second type are those 
from spectra binned up firstly to gather at least 20 counts 
per channel, whereas continuum is modelled again with the 
C-statistic. The third is the case of binned up spectra and 
the model fitted with x 2 statistic. 

The results of the tests are presented in Figs. 6, 7, 8, 
the reference ratios are shown with a solid line, whereas the 
reference average for output bin width equal to 0.1 keV is 
plotted with dots. Tests were performed for different widths 
of output bins, here we present only results for 0.1 keV bins 
since this value corresponds to the ASCA SIS resolution for 
the iron line energy. 

In order to check how the fitted continuum models re- 
produce the reference model, a weighted average of power 
law parameters was calculated for each data set. These re- 
sults are presented in Table 2. Since there is still a small 
tail of the disc line component below 4.5 keV, the power law 
fitted to the reference data (simulated with long exposure) 
is slightly less steep (r = 1.798) than the initial power law 
used in simulations (r = 1.8). All models fitted with the 
C-statistic give almost correct results, whereas fitting with 
the \ 2 statistic leads to a clearly biased result. A similar ef- 
fect was already studied by the Chandra X-ray Center staff 
(Freeman 2001). 

Non-uniform binning up of data before fitting changes 
the relative influence of different parts of the spectrum on 
the fitted model. For a simple continuum model like power 
law this effect is almost negligible (cf. results from row 2 
and rows 3,4 in Table 2), however, for more complex models 
such a procedure should be applied with some caution. 



2 Regarding to arithmetics, 'binning up' is a special case of pre- 
binning, introduced in Sec. 2.2, done for the ouput bin covering 
exactly the sum of input bins. Nonetheless, there is a clear differ- 
ence in applying these two procedures, the former is used before 
fitting the continuum model, the second uses the results of the 
model fitted to the initial input bins. 
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Table 2. Mean values of the continuum model parameters deter- 
mined for various input data and fitting procedure combinations. 
The first row presents the results obtained for the reference data. 



Binning up 


Statistic 


Normalization 


Index 


none 


both 


23.28±0.01 


1.798±0.001 


none 


C, SHERPA 


23.23±0.18 


1.798±0.006 


biased 


C, SHERPA 


23.12±0.19 


1.794±0.006 


random 


C, SHERPA 


23.26±0.19 


1.800±0.006 


biased 


X 2 , XSPEC 


23.55±0.22 


1.829±0.007 



4.3 Biased vs. unbiased binning up 

The standard, initial binning up method, based on adding 
counts from single channels to get at least a given number 
of counts is biased in this sense that resulting bin widths 
are inversely proportional to the measured flux. Then, on 
average, broader input bins are more frequent for ratios be- 
low 1 than for ratios above 1. In consequence, bin weights 
bi are usually smaller for ratios > 1 than bi for ratios < 1 
and this leads to a biased average shape. To test this effect 
we have prepared 'data/model' data with random binning, 
where binning up pattern for a given spectrum was taken 
from the other spectrum, randomly selected. In this way the 
binning should be non-biased, without correlation between 
bin widths and ratios. These spectra were fitted with C- 
statistic and they are the fourth type of 'data/model' data 
tested with some methods. 



4.4 Inhomogeneous shapes 

Consider a special situation: there are two distinct classes of 
objects, with Fe line average shapes clearly different, and, in 
addition, objects of one class are much brighter than those 
of the second class. Then, if the spectra were taken in similar 
conditions (i.e., with approximately equal exposures), their 
weighted average will be far from the true, physical mean 
for these two classes. We have tested such a case by simu- 
lating 20 spectra with 10 times longer exposures and with 
a reference model different from the basic one, described in 
Sec. 4.1. In this model all parameters are the same as in 
the basic one, only the disc line component is weaker, with 
equivalent width equal to the equivalent width of the Gaus- 
sian component, i.e., 60 eV. Using these data and those from 
basic simulations we have checked how the average depends 
on the relative share of different spectra in entire sample. 



5 DISCUSSION 
5.1 Best method 

The best results are obtained for methods la, lb and Ic ap- 
plied to non-binned data, as presented in the bottom panel 
of Fig. 6. Only these methods almost perfectly reproduce the 
shape of the iron line on both its sides. There is a spread of 
results for higher energies, above 7 keV, but these discrep- 
ancies are symmetric and appear due to small statistics in 
this energy range. Methods la and lb used for data binned 
up also work well, but here both red and blue wings of the 
line are slightly but systematically underestimated, as can 




Figure 6. Averaging methods applied to data/model ratios ob- 
tained for non-binned spectra and continuum model fitted with 
C-statistic. Methods are denoted as in Table 1. 



be seen from Fig. 7. For higher energies the spread of re- 
sults is damped due to initial grouping of the input bins, 
however, owing to the same fact, spurious maxima are pro- 
duced around 9.2 and 9.8 keV. Then, even with the best 
methods, averaging binned up data cannot be considered to 
be reliable for the higher energy boundary of ASCA SIS 
range. The third procedure from this group, method Ic, us- 
ing bin weights bi , fails for binned up spectra; this behaviour 
is explained in Sec. 5.2. 

Methods based on averaging with the standard weights 
Wi (Ila-IIe) obviously cannot work properly for a continuum 
modelled with the C-statistic. This is clearly seen in the up- 
per part of Figs. 6 and 7, where all these procedures fail 
completely to reproduce the continuum slope. For the same 
reason, methods using the probability functions do not work 
well for continuum fitted with the \ 2 statistic, this is illus- 
trated in the lower part of Fig. 8. 

The 'standard' method, He, applied to binned up data 
(Fig. 8), can be used as a crude approximation of the average 
shape. The red and blue wings of iron line are here more 
underestimated than is observed for methods Ia-Ic in Fig. 7. 
Also the line peak is not well reproduced and the distortions 
induced by binning up are present in the high energy end 
of the spectrum. Compared with method lib, which uses 
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Figure 7. Averaging methods applied to data/model ratios ob- 
tained for binned up spectra and continuum model fitted with 
C-statistic. 



Figure 8. Averaging methods applied to data/model ratios ob- 
tained for binned up spectra and continuum model fitted with 
X 2 statistic. Correction applied to the results of method Ha is 
explained in the text, see Sec. 5.1. 



all the input bins overlapping with the output bin (with 
gi correction for boundary overlap), the 'standard' method 
results exhibit a larger spread for higher energies. This is 
simply the consequence of neglecting input bins with centres 
lying outside the output bin. Similarly to method Ic applied 
to binned up data (bottom part of Fig. 7), method lie, using 
bin weights, produces a distorted shape for binned up input 
data (top panel of Figs. 7,8). 

Results of two procedures based on prebinning, Ha and 
III, are shown in the middle panel of Figs. 6-8. Due to prebin- 
ning, method Ha is less sensitive to the improper accuracy 
weighting using Wi weights than all procedures averaging 
directly the input ratios (Ilb-IIe). This effect is obviously 
stronger for single instrumental bins (Fig. 6), but appears 
clearly also for rebinned data (Fig. 7,8). 

The arithmetic averaging with prebinning, method III, 
appears to work quite well for data obtained with the C- 
statistic, with, however, a larger spread of results observed 
for higher energies in the case of non-binned data (middle 
panel of Fig. 6) . This spread illustrates the effect of neglect- 
ing any accuracy weights. The overall agreement between 
arithmetic and weighted averages (the proper ones, using 
the probability function) is understandable, due to the fact 
that the longest and the shortest exposure times used in 



simulations differ only by factor of 2.5, thus the accuracy 
weights for different spectra also do not differ strongly. 

The middle part of Fig. 8 shows that the method with 
prebinning and standard accuracy weights Wi, Ha, behaves 
quite differently from its counterpart without prebinning, 
method lib (upper panel of this Figure). This can be ex- 
plained by two effects. The first is an improper continuum 
model, fitted with \ 2 statistics. Using the results presented 
in rows 3 and 5 of Table 2, we have corrected the average 
data/model ratios with a scaling function equal to the ra- 
tio of these mean power law models. The middle panel of 
Fig. 8 presents the results of this correction compared to 
the results of procedure Ha obtained for a continuum model 
fitted with C-statistic (i.e., the same as in the middle panel 
of Fig. 7). Now these results are in agreement and the rest 
of difference between methods lib and Ha, applied to 'x 2 ' 
input data, comes from the fact that the latter procedure is 
less affected by the imperfect accuracy weighting by widths 
Wi. This effect can be estimated from comparison between 
the results of methods lib and Ha shown in the upper and 
middle panels, respectively, of Fig. 7. 
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Figure 9. Test of averaging for biased (correlated) binning up. 
Upper panel: averaging results as in the lower panel of Fig. 7. 
Middle panel shows the ratios between these results and the ref- 
erence shape. The ratios of arithmetic means of the overlap Si for 
data/model ratios above and below 1 are presented in the bottom 
panel. 



5.2 Biased vs. unbiased binning up 

The average iron line shape for Syl nuclei, presented in 
Lubinski & Zdziarski (2001), was obtained with the pro- 
cedure lid applied to initially binned up spectra. This is the 
standard manner used in X-ray spectral fitting with a \ 
statistic to group single input bins to bins with at least 20 
counts. As stated in Sec. 4.3, such a binning up is biased. 
However, it is hard to see how important this effect is in 
the situation where real data of unknown average shape are 
used. A quantitative estimate of the influence of this bias on 
the results of averaging can now be done, for data simulated 
using a known reference model. 

In the upper panel of Fig. 9 the average spectral shapes 
obtained with methods la, lb and Ic for binned up 'C- 
statistic' data are presented again as in Fig. 7. To show more 
clearly the differences between these methods, the ratios be- 
tween their results and the reference average were plotted in 
the middle part of Fig. 9. The bottom panel of this Figure 
illustrates the ratio of the arithmetic averages of the input 
and output bins overlap, Si, determined for data/model ra- 
tios above 1, 5>i, and below 1, <5<i. The deviation between 



Figure 10. The same test as in Fig. 9 done for random (un- 
corrected) binning up. Please note the change of the artificial 
features above 9 keV, in comparison with the results shown in 
Fig. 9, where different binning up pattern was applied. 



the results of method using bin weights, lb, and the refer- 
ence results evidently follows the departure of <5>i/<5<i from 
unity. Fig. 10 presents the results of the same test applied 
to input data binned up with random binning, i.e., where 
the binning pattern for a given spectrum was taken from 
the other spectrum. Now all tested methods are in concor- 
dance, the procedure using bin weights reproduces the spec- 
tral shape equally well. The <5>i/5<i ratios are close to unity, 
but with a larger spread than the same ratios obtained for 
non-binned data, shown for comparison in the bottom part 
of Fig. 10. 

The results of method lid, i.e., that employed by 
Lubinski & Zdziarski (2001), are presented in the upper 
panel of Fig. 8, The correction introduced in Eq. (14) looks 
significant for higher energies (> 6 keV), where gi values are 
larger. However, as we have tested, the differences between 
the results of methods lie and lid disappear when these pro- 
cedures are applied to the data binned up in non-biased way. 
Therefore, correction (14) manifests itself mainly by scaling 
the ratio errors to larger, proper values, almost not affecting 
the ratios themselves. 

The scale of discrepancies between the results of pro- 
cedures using bin weights (Ic,IIc,IId), applied to the data 
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Figure 11. Test for sample with subsamples corresponding to 
two distinct spectral shapes. Upper panel shows the basic (A) and 
alternative (B) reference models together with averages obtained 
with method la for 100 basic (A) and 20 alternative (B) spec- 
tra. In the lower panel average reference models are presented, 
one with weights proportional to the number of spectra in given 
subsample, and second, with equal weights. The average for the 
entire sample of 120 spectra are performed with methods using 
the accuracy weights (la, circles) and without it (III, dots). 

with biased binning up, and the reference results, shown in 
Figs. 7,8, seems to be large. Nevertheless, the spectra tested 
here are extremely weak, and the average spectral shape pre- 
sented in Lubihski & Zdziarski (2001), obtained for better, 
on average, spectra of Syl nuclei, is only moderately affected 
by this mistake. 



Figure 12. Averages with and without bin weights. Results of a 
given method arc presented as the ratio between them and the 
results of the reference averaging. The upper panel shows the 
results for the output bin width equal to 0.1 keV, in the lower 
panel results for 0.05 keV output bins are presented. For clarity, 
results arc grouped (arithmetically) to 0.2 keV bins. 

dition. Afterwards, the weighted averages of approximately 
homogeneous subsamples can be averaged arithmetically. 

In the case of Syl nuclei the spectra are known to be dif- 
ferent for various objects in various spectral states and usu- 
ally the values of their physical parameters are distributed 
over a wide range. However, there may exist samples of ob- 
jects exhibiting clearly distinct physical character. Then the 
average over the entire sample will not describe any real ob- 
ject but, again, this can be easily checked by comparing the 
results for subsamples. 



5.3 Inhomogeneous shapes 

The results of averaging with method la applied to non- 
homogeneous input data are presented in Fig. 11. Two ref- 
erence models are clearly distinct; the blue wings of the disc 
line component in these models are especially different. As- 
suming that the real composition of the studied objects is 
in proportion 1 to 0.2 to the advantage of the stronger disc 
line model, we should expect the mean to be closer to this 
model. However, we have assumed (see Sec. 4.4) that the ob- 
jects with a weaker disc line component are about 10 times 
brighter than the rest of the sample. Then the weighted aver- 
age (circles in the lower panel of Fig. 11) appears to be much 
closer to the shape of the brighter sources in comparison to 
the expected mean (dotted line in this Figure). This disad- 
vantage can be removed with the use of arithmetic average, 
the results of which are shown with dots in the lower part of 
Fig. 11. As already mentioned in Sec. 5.1, arithmetic aver- 
aging leads to a larger spread in results, thus it is advisable 
to check the homogeneity of entire sample by dividing it to 
subsamples, grouped accidentally or according to some con- 



5.4 Incomplete information 

The main reason of using bin weights 6; instead of prebin- 
ning in Lubihski & Zdziarski (2001) was the aim of having 
the possibility of studying all discrete spectral features in the 
average. There was also a second reason: there were many 
spectra, especially those observed with SIS1 spectrometer, 
which, due to the 'worse quality' flag, were truncated far be- 
low the upper instrumental limit (10 keV), at about 7 or 8 
keV. In this case prebinning applied only to completely cov- 
ered output bins leads to the loss of some information. To 
avoid this, one can treat input bins only partially covering 
the output bin as representative of the whole output bin, but 
now their share to the average is somewhat overestimated. 
On the contrary, usage of bin weights automatically reduces 
the significance of such incomplete information. 

To better study the issue of incomplete information we 
have applied procedures with and without prebinning to the 
data, where half, chosen randomly, of the input bins was ne- 
glected in calculating the average. We have tested the best 
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methods, i.e., those based on the product of probability den- 
sity functions, using the data with random binning up, since 
for the non-binned data the effect of bin weighting is re- 
duced. The results of these tests are shown in Fig. 12, for 
two widths of the output bins, 0.1 and 0.05 keV. In the 
case of broader output bins there is no qualitative differ- 
ence between the results of methods la and lb, without bin 
weights, and method Ic, applying b» weights. For narrower 
output bins, in the region of Fe line peak (at 5.5-7 keV) the 
results of method Ic are more consistent with the reference 
average than the results of method lb, whereas method la 
still does not seem to be clearly worse than method Ic. In 
general, bin weighting is advisable in the situation when one 
is going to study the most discrete features of the spectrum. 
However, for weak spectra the advantage of better resolution 
may not prevail over the benefit of higher accuracy coming 
from prcbinning. 

Application of more advanced bin weighting with ad- 
ditional bi/Aa weights (see Sec. 2.1) is also shown in the 
lower part of Fig. 12, but it does not change the results 
significantly. Therefore, there is no need to apply complex 
bin weighting for data collected with instruments of rather 
limited spectral resolution. 

The simultaneous fitting of many spectra, with the same 
line shape parameters but with different continuum models 
and line normalizations, can be considered as an alternative 
to fitting the spectrum obtained from the average shape. 
However, such an approach is inefficient in the case where 
there are more than a few dozen spectra. Handling a large 
set of input data and a much larger set of parameters associ- 
ated with them is impractical for any fitting software due to 
time consuming computations and problems in assigning ac- 
curacy to the fitted values. On the other hand, the spectrum 
with the average shape can be easily fitted when the model 
incorporates some component describing the instrumental 
resolution, e.g., gsmooth model in XSPEC. 



6 CONCLUSIONS 

We have studied the problem of averaging the spectral 
shapes in the case of weak X-ray spectra. The method ap- 
plied in Lubinski & Zdziarski (2001) was substantially im- 
proved by more correct treatment of the Poisson character 
data and better modelling of the reference continuum. Vari- 
ous alternative approaches used to obtain the average spec- 
tral shape were tested with simulated data. 

The reference average is correctly reproduced only by 
the methods based on the description of the data uncer- 
tainty by probability density functions (Eqs. 17-19). Among 
them, the method applying prebinning, i.e., summing the 
number of counts within the output bin (Eq. 18), is the sim- 
plest. Compared to the prebinning method, the procedure 
incorporating weights associated with the input and output 
bin overlap (Eq. 19) works better only in the case where the 
output bins are narrow and, accordingly, it is recommended 
for data taken with instruments of a very good spectral res- 
olution. 

There is no need to initially bin up the averaged data 
when the methods applying the probability functions are 
used and the continuum model is fitted with the C-statistic. 



Only in this case is all of the information collected with the 
single, narrow channels of the instrument taken into account. 
Although the best methods also work well for spectra binned 
up in a non-biased way, the resulting spectral shape can be 
distorted due to the loss or mixing of information contained 
in single bins. 

The usage of the \ 2 statistic in modelling of weak spec- 
tra leads to biased results and deformed average shape. 

The arithmetic average used for prebinned data quite 
fairly reproduces the spectral shape, and can be recom- 
mended for approximately equally accurate input data or 
for testing the homogeneity of the averaged sample. 
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APPENDIX A: AVERAGE FLUX 

Energy flux is defined as the ratio of the number of photons 
Ni collected in a given channel i and the product of the en- 
ergy range of this channel Ae;, exposure time T and effective 
(for the interesting energy range) area A of the instrument 



(fh 



Ni 



ATAe t 



(Al) 



Neglecting the relatively small errors of T, Aa and A 
the accuracy of the flux is given by 



A</>* 



ATAei 



(A2) 



In the case of no background, for Poisson distributed 
data, ANi is usually approximated by \fTN~i, hence the stan- 
dard weighted average of fluxes measured for n k input bins 
overlapping with the output bin has the form 



(ATAd 



fk = 



5 v <» j 



Ni 



ATAei 



AT^Ae, 



(ATAe^ 



(Ae t ) 2 
Ni 



(A3) 



Again using the equality Ni = (f)iAeiAT and inserting 
the output bin width AE k in the numerator and denomina- 
tor we obtain 



E 



Ae, 

aeT 



E 



Aei 



(A4) 



AE k (f)i 



This is the harmonic mean, weighted with bin weights. 
In the situation where flux is approximately constant within 
the output bin, the factor l/{f)i can be taken out of the 
sum in denominator and we can insert (/), in the sum in 
the numerator. Therefore, for the function of flux type, the 
weighted average given by (11) behaves in some sense simi- 
larly to the bin weighted average, Eq. (2) . 

The harmonic average is commonly recommended for 
the quantities of a relative character, defined as a ratio. Nev- 
ertheless, we cannot test this approach since for weak spectra 
the flux measured for narrow input bins is often equal to 0. 



APPENDIX B: POISSON DISTRIBUTION FOR 
UNKNOWN BACKGROUND 

The X-ray background spectra are often weak and the back- 
ground rate cannot be estimated with good accuracy. This 
uncertainty obviously affects the accuracy of the determined 
net source counts. There have been proposed solutions to 
this problem, modifying the Poisson distribution used to de- 
scribe the source + background data. Below we compare the 
results of two of such methods with the distribution given 
by equation (15). The first procedure, classical (Rolke & 
Lopez 2001), is based on the likelihood ratios, the second, 
Bayesian (Loredo 1992), uses marginalization with respect 
to the background rate. For simplicity's sake, we assume here 



that the source and background counts were taken during 
the same time and for the same size of the signal region. 

Rolke & Lopez (2001) developed their method for set- 
ting confidence intervals for small signals in the presence 
of background. To find the confidence region they use the 
likelihood ratio test statistic A in the form 

max{l(X ,ri; N,B) : rj ^ 0} 

A(A ; x, y) = r - t (Bl) 

max{l(X, V ; N, B) : A > 0, 77 > O} 

where A is the tested null hypothesis value of the net 
source counts A. The likelihood function for A source and r\ 
background counts, given the observed N source counts and 
B background counts, is expressed as the product of two 
Poisson distributions 



l(X,V,N,B) 



N\ 



Bl 



(B2) 



As can be expected, the probability distribution derived 
from the likelihood ratio test is broader than the distribution 
obtained for the case of fixed background (Eq. (15)). This is 
the effect of calculating the likelihood function with variable 
r/, what corresponds to taking the background uncertainty 
into account in this method. 

The shape of the likelihood ratio A(A), normalized to 
obtain maximal value equal to 1, is shown in Fig. Bl. In the 
case, when the measured background is equal to 0, broaden- 
ing affects the distribution only in the lowest, close to zero, 
A region. (For B=0 we have modified the Rolke & Lopez 
formula for r\ maximizing the numerator in Eq. (Bl), using 
the relation r/ max — max(Q,(N — 2Ao)/2).) The mode for 
this distribution is equal to N, as in the case of distribu- 
tion obtained for the known background case. However, the 
mean for A(A) distribution is smaller or larger than N + 1, 
depending on N and B. 

The Bayesian method of inferring the signal strength 
in the situation of imprecisely measured background is pre- 
sented in Loredo (1992). First we quote the formula for the 
posterior probability density derived by Loredo for the case 
when the background counts number r) is known 



p(\;N, V ) = C 



(A + ??) 



N„-(X+v) 



N\ 



(B3) 



where the normalization constant C is equal to 
{ Silo r t e T '/^-} ■ The mode for this distribution is equal 
to N and the mean, A, equals to N + 1. 

For the unknown background case, the nuisance param- 
eter, background rate r/, is eliminated through marginaliza- 
tion, i.e., through finding the integral of Eq. (B3) over rj. The 
posterior probability distribution is calculated using the ex- 
pansion of the binomial (A + rf) N and has the form 



p(\;N,B)=J2 C > 



Xe~ 



with coefficients d given by formula 



Ci = 



ni (N+B-i)\ 
Z (N-i)\ 

~j (N+B-j)T 
2^j=0 A (N-j)\ 



(B4) 



(B5) 



The above formula is interpreted in terms of a weighted 
average of the posterior densities calculated attributing 0, 




Figure Bl. Probability distributions for net source counts A ob- 
tained with the use of different methods. N and B are the mea- 
sured source and background counts, respectively. Solid line shows 
distribution obtained for the case, when the background rate is 
assumed to be known (Eq. (15)). Dashed line presents the dis- 
tribution based on the likelihood ratio test, Eq. (Bl). Baycsian 
probability distribution, derived from the marginalization with 
respect to the background rate, Eq. (B4), is plotted with the dot- 
ted line. 



1, N events to the signal. The C; weights are the proba- 
bilities that i of the events observed on-source are from the 
source, provided that B counts are measured off-source. 

The Bayesian probability distribution, computed from 
equation (B4) and normalized to have maximum equal to 
1, differs clearly from two other distributions shown in Fig. 
Bl. Its right tail is more extended but the major differ- 
ence is the shift of its maximum towards lower values of A. 
Depending on how N compares to B, this shift correlates 
inversely with the measured background counts, exceeding 
one unit for B — and approaching for B — N. Using 
the relation J"°° \ l e~ x d\ — i\ it is easy to show that the 
mean A for the probability density given by (B4) is equal 
to J"^_ d(i + 1)/ y^ H - Q Ci- Hence, taking into account the 
background uncertainty in the Bayesian procedure leads to 
lowering of the mean value, with the decrease depending on 
the measured background counts. 

This paper has been typeset from a TgX/ MJrjX file prepared 
by the author. 



